T    = 1.0d4
y    = 0.0d0
nH   = 1.0d-2
G0   = 1.0d-12

sigma0 = 6.3d-18

RC = hiirca_hui(T)
CI = hici_hui(T)

Rn = (CI+RC)*nH
Qn = -( G0 + (CI+2.0d0*RC)*nH + (CI+RC)*nH*y )
Pn = RC * nH * (1.0d0+y)
dn = sqrt(Qn*Qn - 4.0d0*Rn*Pn)

Ri = -(CI+RC)*nH
Qi = CI*nH - G0 - (CI+RC)*nH*y
Pi = G0 + CI*nH*y
di = sqrt(Qi*Qi - 4.0d0*Ri*Pi)

x0i = 0.5d0
x0n = 1.0d0 - x0i

tcn = 2.0d0 / dn
tci = 2.0d0 / di

xsn = (Qn + 2.0d0*Rn*x0n) / dn
xsi = (Qi + 2.0d0*Ri*x0i) / di




N=10000
t_min = 0.0d0
t_max = 5.0d0

t = t_min + dindgen(N)/(N-1) * (t_max - t_min) 


xn = 1.0d0 / (2.0d0 * Rn) * $
  (-Qn - dn * (tanh(t) - xsn)/(1.0d0 - xsn * tanh(t)) )

xi = 1.0d0 / (2.0d0 * Ri) * $
  (-Qi - di * (tanh(t) - xsi)/(1.0d0 - xsi * tanh(t)) )




altay_plot, [0], [0], xrange=[t_min, t_max], yrange=[-4,0.2], $
  xtitle=TexToIdl('t/t_c')

altay_oplot, t, alog10(xn), color=0

altay_oplot, t, alog10(xi), color=0

altay_oplot, t, alog10( 1.0d0 - xn ), color=250, linestyle=2
altay_oplot, t, alog10( 1.0d0 - xi ), color=250, linestyle=2

end
